2020_02_28_lettuce
library(rethinking)## Loading required package: rstan## Loading required package: StanHeaders## Loading required package: ggplot2## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)## Loading required package: parallel## Loading required package: dagitty## rethinking (Version 1.93)##
## Attaching package: 'rethinking'## The following object is masked from 'package:stats':
##
## rstudentlibrary(tidyverse)## ── Attaching packages ───────────────────────────────── tidyverse 1.3.0 ──## ✓ tibble 2.1.3 ✓ dplyr 0.8.3
## ✓ tidyr 1.0.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.4.0
## ✓ purrr 0.3.3## ── Conflicts ──────────────────────────────────── tidyverse_conflicts() ──
## x tidyr::extract() masks rstan::extract()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x purrr::map() masks rethinking::map()Here is your much-awaited dataset for our upcoming meeting on the 21st of February, 2020. The data came fresh off the greenhouse on 2/8/2020 from a night break experiment. Thirty-two accessions from a lettuce F6 recombinant inbred line population were planted in the same greenhouse under 10hr light/14hr dark diurnal cycles. Three treatments are implemented:
- “Nightbreak”: the plants are grown on a bench surrounded by plastic blackout curtains. These plants receive a one-hour night break treatment at 12am every day (meaning the lights turn on in the middle of the night for an hour) in addition to the 10L/14D lighting.
- “Control”: the plants are grown on a bench surrounded by plastic blackout curtains. 10L/14D lighting.
- “Control_NoCurtain”: the plants are grown on a bench without any curtains. 10L/14D lighting.
The goals of the experiment are: a. to see if night break induces earlier flowering in lettuce; b. if so, do different lettuce genotypes respond to night breaks differently; and c. which one(s) of the five candidate loci is/are associated with differential responses.
How to interpret the phenotype: Phenotype is recorded in the “Score” column. The different scores represent different developmental stages: 1: rosette 2: bolted (elongation of the main stem) 3: budding 4: first flower 5: first mature seed head
Aaaaand finally here are your questions! Q1: a. Load the dataset. Look for column “Score” for the response variable we are interested in. A developmental score of 1 or 2 indicates vegetative growth, while a score of 3, 4, or 5 indicates reproductive growth. Create a “Reproduction” column with values 0 and 1, where 0 indicates vegetative growth and 1 indicates reproductive growth.
data <- read_csv("Nightbreak_02_08_20_Rclub.csv")## Parsed with column specification:
## cols(
## RIL = col_character(),
## Treatment = col_character(),
## Rep = col_double(),
## Plot = col_double(),
## Date = col_character(),
## loc1 = col_character(),
## loc2 = col_character(),
## loc3 = col_character(),
## loc4 = col_character(),
## loc5 = col_character(),
## Score = col_double()
## )head(data)## # A tibble: 6 x 11
## RIL Treatment Rep Plot Date loc1 loc2 loc3 loc4 loc5 Score
## <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 100 Control 1 1 2/8/20 P P A P A 2
## 2 2 Control 1 2 2/8/20 A P P A A 2
## 3 137 Control 1 3 2/8/20 P A P P A 1
## 4 172 Control 1 4 2/8/20 P P P A A 2
## 5 188 Control 1 5 2/8/20 P A P P P 2
## 6 PI Control 1 6 2/8/20 P P P P P 3table(data$RIL, data$Treatment)##
## Control Control_NoCurtain NightBreak
## 100 2 2 2
## 105 2 2 2
## 108 2 2 2
## 112 2 2 2
## 114 2 2 2
## 137 2 2 2
## 14 2 2 2
## 146 2 2 2
## 154 2 2 2
## 161 2 2 2
## 169 2 2 2
## 172 2 2 2
## 175 2 2 2
## 179 2 2 2
## 188 2 1 1
## 19 2 2 2
## 2 2 2 2
## 201 2 2 2
## 21 2 2 2
## 222 2 2 2
## 228 2 2 2
## 229 2 2 2
## 235 2 2 2
## 24 2 2 2
## 25 2 1 1
## 33 2 2 2
## 53 2 2 2
## 55 2 2 2
## 58 2 2 2
## 66 2 2 2
## Arm 2 2 2
## PI 2 2 2data <- data %>%
mutate(reproduction=ifelse(Score>2, 1, 0))- (optional) Take a look at columns “loc1” through “loc5”. The levels “A” or “P” indicate which parental allele the RIL has inherited. Can you think of a reason why there are 32 lines in this experiment?
2^5## [1] 32Q2: Using the “Reproduction” column you created in Q1a as the response variable, fit a simple model with effects of genotype (“RIL”) and treatment (“Treatment”) on the rate of transitioning to reproductive growth by 2/8/20. (Things you might want to consider: appropriate likelihood function, intersection term).
First fit a model with Treatment only
data <- data %>%
mutate(RIL_i = as.numeric(as.factor(data$RIL)),
NoCurtain = ifelse(Treatment=="Control_NoCurtain", 1L, 0L),
NightBreak = ifelse(Treatment=="NightBreak", 1L, 0L))
dat2a <- data %>% select(reproduction, NoCurtain, NightBreak)
fm2a <- ulam(alist(reproduction ~ dbinom(1, p),
logit(p) <- a + b_nc*NoCurtain + b_nb*NightBreak,
a ~ dnorm(0, 1.5),
b_nc ~ dnorm(0,1),
b_nb ~ dnorm(0,1)),
data=dat2a,
cores = 4,
chains = 4,
log_lik=TRUE)precis(fm2a)## mean sd 5.5% 94.5% n_eff Rhat
## a -1.3684195 0.2809390 -1.8362061 -0.9395667 607.5517 1.007448
## b_nc -0.2799247 0.4174428 -0.9235534 0.4000005 674.5954 1.003021
## b_nb 1.9834379 0.3662811 1.4359085 2.5824523 679.5032 1.005163pairs(fm2a)
traceplot(fm2a)
trankplot(fm2a)
coef(fm2a)## a b_nc b_nb
## -1.3684195 -0.2799247 1.9834379inv_logit(coef(fm2a)["a"])## a
## 0.2028753inv_logit(coef(fm2a)["a"] + coef(fm2a)["b_nb"])## a
## 0.6490847Night breaks promote reproduction probability
Now try it with RIL. Probably this should be donw pooling info across RILs but not doing that now.
dat2b <- data %>% select(reproduction, RIL_i, NoCurtain, NightBreak)
fm2b <- ulam(alist(reproduction ~ dbinom(1, p),
logit(p) <- a[RIL_i] +
b_nc*NoCurtain +
b_nb*NightBreak +
i_nc[RIL_i]*NoCurtain +
i_nb[RIL_i]*NightBreak,
a[RIL_i] ~ dnorm(0, 1.5),
b_nc ~ dnorm(0, 1),
b_nb ~ dnorm(0, 1),
i_nc[RIL_i] ~ dnorm(0, .5),
i_nb[RIL_i] ~ dnorm(0, .5)),
data=dat2b,
cores = 4,
chains = 4,
log_lik=TRUE)traceplot(fm2b, as=FALSE)## Waiting to draw page 2 of 7
## Waiting to draw page 3 of 7
## Waiting to draw page 4 of 7
## Waiting to draw page 5 of 7
## Waiting to draw page 6 of 7
## Waiting to draw page 7 of 7
precis(fm2b, depth=2)## mean sd 5.5% 94.5% n_eff Rhat
## a[1] -1.513511449 0.9053477 -2.9574742 -0.10964049 3415.186 0.9992541
## a[2] -1.531490572 0.9096282 -3.0393241 -0.09293946 2766.393 0.9990140
## a[3] -2.328239036 1.0218807 -3.9627221 -0.77726757 2172.867 0.9988177
## a[4] -0.786342529 0.8846455 -2.2461891 0.58659280 3060.438 0.9987300
## a[5] -2.314416408 0.9906888 -3.9913930 -0.78990458 2955.146 0.9990294
## a[6] -2.329416752 0.9881313 -3.9118923 -0.83385316 3149.466 0.9999013
## a[7] 2.201275809 0.9671207 0.7560101 3.79382037 2952.022 1.0000600
## a[8] 2.205695208 0.9871531 0.6472986 3.85478246 3095.989 0.9989667
## a[9] -0.803567856 0.8873465 -2.1988309 0.60451484 4091.916 0.9984762
## a[10] 2.159343184 0.9303570 0.7296120 3.64486085 2805.630 0.9996754
## a[11] -0.778614950 0.8886038 -2.2298829 0.62816630 2923.399 0.9994641
## a[12] -0.807033541 0.8750642 -2.1848168 0.58990148 2612.720 0.9987794
## a[13] -0.825370948 0.8937957 -2.2464411 0.59078410 2508.108 0.9990964
## a[14] -0.800126880 0.8568508 -2.2229429 0.54448881 3821.607 0.9981235
## a[15] -0.045965414 0.9727623 -1.6227124 1.50971440 3340.204 0.9990684
## a[16] -0.787698969 0.8799377 -2.2116559 0.56357827 2356.167 0.9996631
## a[17] -0.787199172 0.8607496 -2.1561030 0.55454757 2926.142 0.9992964
## a[18] -0.766607280 0.8842022 -2.2278730 0.69292667 2917.478 1.0013686
## a[19] -2.365884408 0.9622792 -3.9721218 -0.88546106 3467.740 0.9989408
## a[20] -1.515210447 0.9021454 -2.9634031 -0.09589020 2672.041 0.9985085
## a[21] -0.799482077 0.8940985 -2.2002786 0.61357345 2855.338 0.9986670
## a[22] -2.306206799 0.9820828 -3.9398599 -0.75087240 2402.785 0.9994058
## a[23] 2.201332871 1.0154349 0.6769680 3.87670461 3841.979 0.9992756
## a[24] -0.806807346 0.8821051 -2.2165338 0.58877544 3477.490 0.9993531
## a[25] -0.943549209 0.9235160 -2.4237154 0.55600723 2669.832 0.9997891
## a[26] -2.311572487 1.0289183 -4.0531471 -0.79125455 2653.576 1.0001946
## a[27] -0.794526866 0.8998186 -2.2555147 0.59879446 2907.020 0.9984286
## a[28] -1.516998264 0.9090782 -3.0296035 -0.07855433 2418.267 1.0014764
## a[29] -2.336684963 0.9779970 -3.9513996 -0.84028024 2714.666 0.9993002
## a[30] -0.785537610 0.8904814 -2.2766513 0.62214688 2798.109 1.0004535
## a[31] -2.365032519 1.0092728 -4.0364525 -0.86260553 2846.423 1.0003047
## a[32] 2.203717128 0.9928840 0.6453905 3.83422473 3221.755 0.9990227
## b_nc -1.461095658 0.4473172 -2.2138761 -0.75628451 2356.909 0.9984194
## b_nb 1.799079807 0.3631064 1.2156199 2.37878936 1706.030 1.0025933
## i_nc[1] -0.030388482 0.4746775 -0.7988804 0.73573395 3142.079 0.9994152
## i_nc[2] -0.036604057 0.4845529 -0.8212190 0.70768886 4460.628 0.9992232
## i_nc[3] -0.016019849 0.4995299 -0.8321545 0.78814566 4059.139 0.9990353
## i_nc[4] -0.066080602 0.4894878 -0.8476785 0.69956371 2515.776 0.9990334
## i_nc[5] -0.019836458 0.5044797 -0.8224641 0.78861174 3655.672 0.9991428
## i_nc[6] -0.008471087 0.4923143 -0.8023173 0.78824561 2672.513 0.9990305
## i_nc[7] 0.166824028 0.5075839 -0.6344034 0.96195718 3629.027 0.9990174
## i_nc[8] 0.157021952 0.4878462 -0.6117654 0.94676182 3375.007 0.9998238
## i_nc[9] -0.062249139 0.4826179 -0.8368033 0.70597007 2477.363 0.9986126
## i_nc[10] 0.150809628 0.4914341 -0.6590687 0.93449341 3478.837 0.9986285
## i_nc[11] -0.061566214 0.5076639 -0.8602698 0.74330480 2974.488 0.9986786
## i_nc[12] -0.068733771 0.4978988 -0.8619567 0.74148771 2845.342 0.9986862
## i_nc[13] -0.070907597 0.4889094 -0.8639211 0.71227173 3205.711 0.9993335
## i_nc[14] -0.058162518 0.5169410 -0.8664763 0.77376184 3263.033 0.9992405
## i_nc[15] -0.052647159 0.5021671 -0.8480280 0.75345010 2972.775 0.9986718
## i_nc[16] -0.060113157 0.5129813 -0.9033912 0.73155620 4965.434 0.9986576
## i_nc[17] -0.061967188 0.4982477 -0.8384984 0.69742916 3174.896 0.9990903
## i_nc[18] -0.064601816 0.4753955 -0.8027760 0.67777649 3028.033 0.9983097
## i_nc[19] -0.021465629 0.4853945 -0.7946802 0.76718013 2722.543 0.9999162
## i_nc[20] -0.034476559 0.4893756 -0.7942391 0.74682488 3415.669 0.9997953
## i_nc[21] -0.042842433 0.4801999 -0.8061453 0.71815898 3110.034 0.9986689
## i_nc[22] -0.020282222 0.4998211 -0.8349103 0.75682534 3232.972 0.9996115
## i_nc[23] 0.163254583 0.4751780 -0.5901927 0.92460240 2891.293 0.9994749
## i_nc[24] -0.071562327 0.4902698 -0.8426906 0.73126364 3731.281 0.9986923
## i_nc[25] -0.032224630 0.4993519 -0.8107605 0.76309232 3072.595 0.9984978
## i_nc[26] -0.028286124 0.5117882 -0.8353263 0.78786794 3558.199 0.9989224
## i_nc[27] -0.056234459 0.4771904 -0.8210768 0.71461932 3299.372 0.9987108
## i_nc[28] -0.051002121 0.5115588 -0.8572915 0.78692026 3817.466 0.9998203
## i_nc[29] -0.012640750 0.5126681 -0.8571728 0.79765597 4726.750 0.9992283
## i_nc[30] -0.068959462 0.4790384 -0.8447502 0.70231895 2619.255 1.0008454
## i_nc[31] -0.015330912 0.4819963 -0.7750983 0.74486117 3954.618 0.9989296
## i_nc[32] 0.163793035 0.4680138 -0.5756983 0.90287120 3557.668 0.9984449
## i_nb[1] -0.019333710 0.4889802 -0.8215391 0.72366721 3636.266 1.0007879
## i_nb[2] -0.032667914 0.4696069 -0.7731813 0.70115006 2791.465 0.9992627
## i_nb[3] -0.167052037 0.4992810 -0.9746051 0.62978460 2857.589 0.9982153
## i_nb[4] 0.144214095 0.4889962 -0.6895076 0.92196835 2517.590 1.0002884
## i_nb[5] -0.182035837 0.4521319 -0.9064093 0.53368006 2745.134 0.9989976
## i_nb[6] -0.178218039 0.4742781 -0.9447832 0.55477722 3474.378 0.9986948
## i_nb[7] 0.014953622 0.4816249 -0.7359239 0.79785225 3753.746 0.9988961
## i_nb[8] 0.020810560 0.4933788 -0.7862159 0.82052178 2756.127 0.9998826
## i_nb[9] 0.141120649 0.4724785 -0.5904523 0.90355759 3486.947 0.9992673
## i_nb[10] 0.025798362 0.4881785 -0.7355953 0.82256976 3587.553 0.9986180
## i_nb[11] 0.128544232 0.4951204 -0.6779098 0.90076542 3151.730 0.9994870
## i_nb[12] 0.145148832 0.5213247 -0.6934230 1.00024227 4972.545 0.9983697
## i_nb[13] 0.143985022 0.4818559 -0.6277442 0.92591942 2876.482 0.9997406
## i_nb[14] 0.142117612 0.4943731 -0.6734830 0.90910712 3271.333 0.9987387
## i_nb[15] 0.048223519 0.4805785 -0.7060056 0.79508991 3867.197 0.9988779
## i_nb[16] 0.135229146 0.4910003 -0.6589576 0.92562032 3304.149 0.9991042
## i_nb[17] 0.131152978 0.4879140 -0.6292739 0.90561689 3622.675 0.9988634
## i_nb[18] 0.132252996 0.4613578 -0.6166871 0.88159440 3973.873 0.9990162
## i_nb[19] -0.187153676 0.4784037 -0.9624071 0.59158162 3302.967 0.9998225
## i_nb[20] -0.027038186 0.4722429 -0.7638541 0.74350136 3989.402 0.9993926
## i_nb[21] 0.150629118 0.4896246 -0.6451774 0.95178746 3457.173 0.9990643
## i_nb[22] -0.192550850 0.4807391 -0.9793850 0.58555045 3596.781 0.9991299
## i_nb[23] 0.015813601 0.4748413 -0.7210722 0.76887061 3091.648 0.9994998
## i_nb[24] 0.138618157 0.4733826 -0.6484403 0.88467280 2742.717 1.0010410
## i_nb[25] 0.083104216 0.5037001 -0.7161947 0.87063618 3359.241 0.9991885
## i_nb[26] -0.183307926 0.5134516 -0.9841362 0.65501731 4172.173 0.9993981
## i_nb[27] 0.133966010 0.4844314 -0.6493674 0.91503726 3649.225 0.9988449
## i_nb[28] -0.027622467 0.4882822 -0.8058887 0.76299970 3533.817 0.9994261
## i_nb[29] -0.175055339 0.4961503 -0.9724825 0.62537558 3569.176 0.9993246
## i_nb[30] 0.138814689 0.4915381 -0.6348402 0.94704551 3215.567 0.9992382
## i_nb[31] -0.176839930 0.4936670 -0.9735989 0.59763975 3725.328 0.9988709
## i_nb[32] 0.026540425 0.4963812 -0.7597550 0.82028595 3669.548 0.9986951
precis(fm2b, depth=2) %>%
rownames_to_column("parameter") %>%
mutate(parameter_type=str_sub(parameter, 1, 1)) %>%
ggplot(aes(x=parameter, y=mean, ymax=`94.5%`, ymin=`5.5%`)) +
geom_errorbar() +
geom_point() +
facet_wrap(facets=~parameter_type, scales="free_x",ncol=1) +
theme(axis.text.x = element_text(angle=90, hjust=1, size = 8))
Different intercepts for the different RILs, but no evidence of interaction.
Q3: Because we are more interested in the effects of individual loci than the performance of specific genotypes, fit a model with additive effects of the five loci and effect of treatment.
dat3 <- data %>% select(reproduction, NoCurtain, NightBreak, starts_with("loc")) %>%
mutate_at(vars(starts_with("loc")), ~ ifelse(.=="P", 1L, 0L))
dat3## # A tibble: 188 x 8
## reproduction NoCurtain NightBreak loc1 loc2 loc3 loc4 loc5
## <dbl> <int> <int> <int> <int> <int> <int> <int>
## 1 0 0 0 1 1 0 1 0
## 2 0 0 0 0 1 1 0 0
## 3 0 0 0 1 0 1 1 0
## 4 0 0 0 1 1 1 0 0
## 5 0 0 0 1 0 1 1 1
## 6 1 0 0 1 1 1 1 1
## 7 0 0 0 0 0 0 1 1
## 8 0 0 0 1 1 0 0 0
## 9 1 0 0 0 1 1 0 1
## 10 0 0 0 1 0 0 1 1
## # … with 178 more rowsfm3 <- ulam(alist(reproduction ~ dbinom(1, p),
logit(p) <- a +
b_nc*NoCurtain +
b_nb*NightBreak +
b1*loc1 +
b2*loc2 +
b3*loc3 +
b4*loc4 +
b5*loc5,
a ~ dnorm(0, 1.5),
b_nc ~ dnorm(0, 1),
b_nb ~ dnorm(0, 1),
c(b1, b2, b3, b4, b5) ~ dnorm(0,1)),
data=dat3,
cores = 4,
chains = 4,
log_lik=TRUE)traceplot(fm3)
trankplot(fm3)
pairs(fm3)
precis(fm3)## mean sd 5.5% 94.5% n_eff Rhat
## a -3.3159824 0.5342490 -4.1645083 -2.4959511 1207.117 0.9991905
## b_nc -0.4508491 0.4501747 -1.1826645 0.2435275 1868.081 0.9997200
## b_nb 2.2822734 0.4159139 1.6198527 2.9565003 1701.503 0.9984651
## b5 0.9734930 0.3600145 0.4100787 1.5525136 1776.041 0.9995795
## b4 0.2486990 0.3601993 -0.2937959 0.8405605 1830.149 1.0019061
## b3 1.1465042 0.3639615 0.5607967 1.7229856 1714.689 0.9983881
## b2 1.3472405 0.3501848 0.7882069 1.9125313 1791.593 0.9995878
## b1 -0.1771559 0.3468650 -0.7287303 0.3767315 1920.996 0.9990492plot(precis(fm3))
Q4: Now let’s look at some interaction terms. Can you fit a model that takes into account interaction effects between treatment and allele types at the five loci? How do you interpret the output? (I built a somewhat “ugly” model for this question. I’m excited to see what y’all’s models look like.)
fm4 <- ulam(alist(reproduction ~ dbinom(1, p),
logit(p) <- a +
b_nc*NoCurtain +
b_nb*NightBreak +
b1*loc1 +
b2*loc2 +
b3*loc3 +
b4*loc4 +
b5*loc5 +
i_nc1*NoCurtain*loc1 +
i_nc2*NoCurtain*loc2 +
i_nc3*NoCurtain*loc3 +
i_nc4*NoCurtain*loc4 +
i_nc5*NoCurtain*loc5 +
i_nb1*NightBreak*loc1 +
i_nb2*NightBreak*loc2 +
i_nb3*NightBreak*loc3 +
i_nb4*NightBreak*loc4 +
i_nb5*NightBreak*loc5 ,
a ~ dnorm(0, 1.5),
b_nc ~ dnorm(0, 1),
b_nb ~ dnorm(0, 1),
c(b1, b2, b3, b4, b5) ~ dnorm(0,1),
c(i_nc1, i_nc2, i_nc3, i_nc4, i_nc5) ~ dnorm(0,1),
c(i_nb1, i_nb2, i_nb3, i_nb4, i_nb5) ~ dnorm(0,1)
),
data=dat3,
cores = 4,
chains = 4,
log_lik=TRUE)traceplot(fm4, ask=FALSE)## Waiting to draw page 2 of 2
precis(fm4)## mean sd 5.5% 94.5% n_eff Rhat
## a -3.20206160 0.5738069 -4.1567922 -2.30902682 1440.342 1.0011288
## b_nc -1.13186491 0.7309265 -2.2724402 0.02684157 1928.929 0.9999408
## b_nb 1.86162350 0.6442929 0.8296973 2.89968701 1510.263 1.0037985
## b5 0.90847570 0.4588217 0.1722940 1.64847100 1824.613 0.9999606
## b4 0.19550106 0.4660985 -0.5640710 0.95615009 1831.208 1.0038980
## b3 0.99087955 0.4568493 0.2612145 1.71951886 1579.873 0.9999571
## b2 1.30512184 0.4624019 0.5721861 2.07501096 1795.455 0.9990708
## b1 -0.36202761 0.4443506 -1.0614068 0.35388762 1740.088 0.9994106
## i_nc5 0.32568212 0.6316553 -0.6783409 1.35950807 2167.702 0.9986779
## i_nc4 0.04120913 0.6676923 -1.0002293 1.12584010 2398.532 1.0004317
## i_nc3 0.29742951 0.6462575 -0.7426241 1.34074642 1984.749 0.9986551
## i_nc2 0.81750985 0.6754847 -0.2639218 1.91081049 2184.731 0.9988051
## i_nc1 -0.34465373 0.6910576 -1.4350683 0.73902022 2619.572 0.9997956
## i_nb5 0.12110959 0.6061301 -0.8496245 1.07471432 1919.323 1.0001280
## i_nb4 0.12897554 0.6165936 -0.8738151 1.11489907 1964.536 1.0038060
## i_nb3 0.47304185 0.6309728 -0.5336297 1.46997981 1996.184 0.9997948
## i_nb2 -0.13227369 0.6138426 -1.1314829 0.85586888 2458.062 1.0006682
## i_nb1 0.69192460 0.6019424 -0.2730246 1.66142048 2149.760 0.9999572plot(precis(fm4))
Q5: By simplifying the developmental score phenotype into a binary variable that indicates whether a plant has entered reproductive growth, we run the risk of losing potentially important information. Re-fit your favorite model from Q4 with the ordered categorical outcome variable of “Score.” Do you observe any changes in your results? If so, why do you think it happened?
dat5 <- data %>% select(Score, NoCurtain, NightBreak, starts_with("loc")) %>%
mutate_at(vars(starts_with("loc")), ~ ifelse(.=="P", 1L, 0L))
dat5## # A tibble: 188 x 8
## Score NoCurtain NightBreak loc1 loc2 loc3 loc4 loc5
## <dbl> <int> <int> <int> <int> <int> <int> <int>
## 1 2 0 0 1 1 0 1 0
## 2 2 0 0 0 1 1 0 0
## 3 1 0 0 1 0 1 1 0
## 4 2 0 0 1 1 1 0 0
## 5 2 0 0 1 0 1 1 1
## 6 3 0 0 1 1 1 1 1
## 7 1 0 0 0 0 0 1 1
## 8 2 0 0 1 1 0 0 0
## 9 3 0 0 0 1 1 0 1
## 10 1 0 0 1 0 0 1 1
## # … with 178 more rowsfm5 <- ulam(alist(Score ~ dordlogit(phi, cutpoints),
logit(phi) <-
b_nc*NoCurtain +
b_nb*NightBreak +
b1*loc1 +
b2*loc2 +
b3*loc3 +
b4*loc4 +
b5*loc5 +
i_nc1*NoCurtain*loc1 +
i_nc2*NoCurtain*loc2 +
i_nc3*NoCurtain*loc3 +
i_nc4*NoCurtain*loc4 +
i_nc5*NoCurtain*loc5 +
i_nb1*NightBreak*loc1 +
i_nb2*NightBreak*loc2 +
i_nb3*NightBreak*loc3 +
i_nb4*NightBreak*loc4 +
i_nb5*NightBreak*loc5 ,
b_nc ~ dnorm(0, 1),
b_nb ~ dnorm(0, 1),
c(b1, b2, b3, b4, b5) ~ dnorm(0,1),
c(i_nc1, i_nc2, i_nc3, i_nc4, i_nc5) ~ dnorm(0,1),
c(i_nb1, i_nb2, i_nb3, i_nb4, i_nb5) ~ dnorm(0,1),
cutpoints ~ dnorm(0, 1.5)
),
data=dat5,
cores = 4,
chains = 4,
log_lik=TRUE)traceplot(fm5, ask=FALSE)## Waiting to draw page 2 of 2
trankplot(fm5, ask=FALSE)
## Waiting to draw page 2 of 2
precis(fm5)## 3 vector or matrix parameters hidden. Use depth=2 to show them.## mean sd 5.5% 94.5% n_eff Rhat
## b_nc -1.36461870 0.8129094 -2.66796937 -0.02584995 2495.925 1.0000141
## b_nb 1.31235615 0.8499278 -0.07734712 2.65725972 3166.144 0.9992015
## b5 0.18469567 0.8010326 -1.09751294 1.44943899 2397.151 0.9992830
## b4 -0.32922449 0.8027642 -1.56741531 1.00838988 2584.827 0.9997975
## b3 0.23012416 0.7977883 -0.97898000 1.52957173 2055.444 1.0007466
## b2 0.53442484 0.8843423 -0.86369084 1.97109490 2527.667 0.9998631
## b1 -0.12002411 0.8308789 -1.37015042 1.25020822 3362.073 0.9989224
## i_nc5 -0.11024793 0.8893274 -1.55475156 1.29715673 3032.845 1.0008822
## i_nc4 -0.24092744 0.9359842 -1.80776044 1.27162330 3090.932 0.9990354
## i_nc3 0.03364703 0.9688201 -1.49997410 1.54912700 2820.165 0.9993861
## i_nc2 0.14128244 0.9811925 -1.40131002 1.76925999 2672.643 1.0025950
## i_nc1 -0.02878890 0.9120195 -1.44865318 1.44108208 3147.552 0.9984879
## i_nb5 0.61401854 0.8925083 -0.78630883 2.00073861 3715.482 0.9999595
## i_nb4 0.63803136 0.8737669 -0.74570706 1.99913716 3257.582 0.9992951
## i_nb3 0.63323963 0.8855830 -0.80001976 2.00581768 2916.897 1.0012961
## i_nb2 0.53799210 0.9641118 -0.98001782 2.10110636 3479.199 1.0011805
## i_nb1 0.65972872 0.8789137 -0.69748381 2.04565876 3035.606 0.9984355plot(precis(fm5))## 3 vector or matrix parameters hidden. Use depth=2 to show them.
Now there are no effect of the loci. Loci are specific to reproduction rather than overall development?
Q6: Each “Plot” # correspond to a specific spot on a bench. In other words, the same plot # indicates equivalent locations on their respective benches even across different treatments and replicates. Update your favorite model from Q4 or Q5 using hierarchical modeling that allow partial pooling across plots. Compare the models. What do they say and which model do you prefer?
dat6 <- data %>% select(reproduction, NoCurtain, NightBreak, starts_with("loc"), Plot) %>%
mutate_at(vars(starts_with("loc")), ~ ifelse(.=="P", 1L, 0L))
sort(unique(dat6$Plot))## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32dat6## # A tibble: 188 x 9
## reproduction NoCurtain NightBreak loc1 loc2 loc3 loc4 loc5 Plot
## <dbl> <int> <int> <int> <int> <int> <int> <int> <dbl>
## 1 0 0 0 1 1 0 1 0 1
## 2 0 0 0 0 1 1 0 0 2
## 3 0 0 0 1 0 1 1 0 3
## 4 0 0 0 1 1 1 0 0 4
## 5 0 0 0 1 0 1 1 1 5
## 6 1 0 0 1 1 1 1 1 6
## 7 0 0 0 0 0 0 1 1 7
## 8 0 0 0 1 1 0 0 0 8
## 9 1 0 0 0 1 1 0 1 9
## 10 0 0 0 1 0 0 1 1 10
## # … with 178 more rowsfm6 <- ulam(alist(reproduction ~ dbinom(1, p),
logit(p) <- a_bar + z[Plot]*sigma_a +
b_nc*NoCurtain +
b_nb*NightBreak +
b1*loc1 +
b2*loc2 +
b3*loc3 +
b4*loc4 +
b5*loc5 +
i_nc1*NoCurtain*loc1 +
i_nc2*NoCurtain*loc2 +
i_nc3*NoCurtain*loc3 +
i_nc4*NoCurtain*loc4 +
i_nc5*NoCurtain*loc5 +
i_nb1*NightBreak*loc1 +
i_nb2*NightBreak*loc2 +
i_nb3*NightBreak*loc3 +
i_nb4*NightBreak*loc4 +
i_nb5*NightBreak*loc5 ,
z[Plot] ~ dnorm(0, 1),
a_bar ~ dnorm(0, 1.5),
sigma_a ~ dexp(1),
b_nc ~ dnorm(0, 1),
b_nb ~ dnorm(0, 1),
c(b1, b2, b3, b4, b5) ~ dnorm(0,1),
c(i_nc1, i_nc2, i_nc3, i_nc4, i_nc5) ~ dnorm(0,1),
c(i_nb1, i_nb2, i_nb3, i_nb4, i_nb5) ~ dnorm(0,1)
),
data=dat6,
cores = 4,
chains = 4,
log_lik=TRUE)traceplot(fm6, ask=FALSE)## Waiting to draw page 2 of 4
## Waiting to draw page 3 of 4
## Waiting to draw page 4 of 4
precis(fm6)## 32 vector or matrix parameters hidden. Use depth=2 to show them.## mean sd 5.5% 94.5% n_eff Rhat
## a_bar -3.22553855 0.5896173 -4.18513147 -2.28836963 2316.4285 0.9984087
## sigma_a 0.38577990 0.3033213 0.02467026 0.95410480 956.1048 1.0003018
## b_nc -1.15021044 0.6905816 -2.29169987 -0.03358274 2571.0776 0.9998338
## b_nb 1.88120553 0.6329843 0.88468289 2.87650839 2321.7827 0.9984949
## b5 0.91309885 0.4784659 0.17977721 1.71354018 2163.2681 0.9996737
## b4 0.18428139 0.4762249 -0.56983685 0.94755433 2790.5947 0.9987916
## b3 1.00768002 0.4777706 0.24372598 1.77391019 2470.8216 0.9985788
## b2 1.32235272 0.4805014 0.57713580 2.13698333 2481.9108 0.9985999
## b1 -0.42174529 0.4740164 -1.21390515 0.34283798 2338.2097 1.0004988
## i_nc5 0.33343033 0.6942860 -0.79617872 1.42898905 3120.4322 0.9998352
## i_nc4 0.04423745 0.6455524 -0.98218006 1.06514654 3092.4307 0.9993721
## i_nc3 0.30475783 0.6763374 -0.73119744 1.38797602 3048.3032 0.9986035
## i_nc2 0.78872132 0.6699583 -0.28858556 1.84791081 3590.5671 0.9987599
## i_nc1 -0.30409907 0.6857999 -1.41496435 0.77914038 3577.9156 0.9990364
## i_nb5 0.11442120 0.6071274 -0.85961686 1.07800734 2653.7465 0.9986733
## i_nb4 0.16337013 0.6316318 -0.84048736 1.17161341 3174.1875 0.9983739
## i_nb3 0.54704199 0.6036089 -0.42835394 1.50796303 2941.8996 0.9992783
## i_nb2 -0.12372440 0.6227332 -1.10401526 0.89147085 2946.6912 0.9988911
## i_nb1 0.74884570 0.6155998 -0.21847710 1.72967241 2465.1401 0.9994735plot(precis(fm6))## 32 vector or matrix parameters hidden. Use depth=2 to show them.
compare(fm4, fm6)## WAIC SE dWAIC dSE pWAIC weight
## fm4 168.8251 12.65354 0.000000 NA 10.73038 0.6418788
## fm6 169.9922 12.79354 1.167056 1.298259 14.45918 0.3581212Q7 (optional): a. What can we conclude regarding treatment effect?
Night breaks increases reproduction probability. For example, when all loci are “A”:
inv_logit(coef(fm4)["a"])## a
## 0.03908821inv_logit(coef(fm4)["a"] + coef(fm4)["b_nb"])## a
## 0.207438- What can we conclude regarding differential response to nightbreak treatment?
Although some loci affect reproduction, there is no interaction between loci and night breaks.
- If we concluded that there are differential responses to nightbreak across genotypes, which genetic loci contributed to the differential responses, and which ones did not?